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Fragmentation of a sheet by propagating, branching and merging cracks 
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We consider a model of fragmentation of sheet by cracks that move with a velocity in preferred 
direction, but undergo random transverse displacements as they move. There is a non-zero proba¬ 
bility of crack-splitting, and the split cracks move independently. If two cracks meet, they merge, 
and move as a single crack. In the steady state, there is non-zero density of cracks, and the sheet 
left behind by the moving cracks is broken into a large number of fragments of different sizes. The 
evolution operator for this model reduces to the Hamiltonian of quantum XY spin chain, which is 
exactly integrable. This allows us to determine the steady state, and also the distribution of sizes 
of fragments. 

PACS numbers: 75.10.Jm 


I. INTRODUCTION 

There has been huge body of work in engineering and physics literature dealing with the distribution of sizes 
of fragments, when a large piece of solid is broken into much smaller fragments, driven by its applications in many 
industrial processes like mining and milling, in geophysics and environmental science [IHU- The theoretical approaches 
to modelling fragmentation have generally focussed on the distribution of sizes of fragments. It is expected that the 
general form of the distribution would be independent of material details, and may be captured in a simple theoretical 
model based on some general physical principles. In one of the first theoretical attempts in this direction, Mott 
and Linfoot [5] modelled the fragmentation as caused by randomly drawn Poisson-distributed horizontal and vertical 
lines. Grady and Kipp studied different variants of the Mott construction, and their effect on the fragment-size 
distribution 1]. Similar ideas were also used by Gilvarry mm to predict the fragment size distribution. These have 
been studied using the chemical reaction rate equations approach [51. Several random fuse-network type models for 
fracture have been developed, where the material is seen as a network of coupled elastic elements, which have random 
threshold for failure ®iin! . There are also molecular-dynamic simulations of assembly of sub-units, interacting by 
say Lennard-Jones interaction, subjected to high-impact projectile, or expansive stress D3HI51- Models assuming a 
fractal structure of defects have also been studied m- 

In recent years, there has been a lot of concern about climate change, leading to increased melting of the polar ice, 
and consequent rise of sea-levels. To make quantitative predictions, one needs to understand the fragmentation of 
the the polar ice cap into much smaller iceberg fragments. This led to recent study of calving glaciers by Astrorn et 
al |l7]. The precise model studied by Astrorn et al is rather complicated, and has to be solved numerically [lHl H]■ 
In this paper, we propose a simplified model of stochastic growth of cracks on a two dimensional sheet, which was 
inspired by these papers. The cracks have a preferred direction of propagation, but have a random transverse velocity. 
There is a finite probability per unit time that the a crack splits into two, which then propagate independently. If two 
cracks meet, they merge, and evolve subsequently as a single crack. The region left behind by the cracks is broken 
into fragments of different sizes. An example of the fragments generated in our model is shown in Fig. [T] 

This model was studied earlier by Ben Avraham et al EDI 122 as a model of particles with diffusion, aggregation 
and birth. They noted that the detailed balance condition is satisfied, and hence the steady state is easily written 
down, and showed that the spacing distribution between consecutive particles satisfies a linear equation, and hence 
can be determined. The evolution operator for this model reduces to the Hamiltonian of a quantum XY spin model, 
which well known to be diagonalizable exactly in terms of free fermions |22H23| . Viewed as a model of fragmentation, 
geometrical questions about fragment-size distribution and correlations seem quite natural, but these have not been 
discussed before. A model rather similar in spirit was discussed earlier by Inaoka and Takayasu ]26l . but their detailed 
model is actually quite different, and can only be solved numerically. Our model has the advantage of been analytically 
tractable, and we are able determine the exact distribution of fragment sizes. We show that in the limit when the 
splitting rate is small, there is an exact mapping between our model A (defined precisely below) and the directed 
abelian sandpile model EH EH!, which brings out the self-organized origin of the divergence in the fragment-size 
distribution for small sizes. In fact, this distribution is the same as that encountered in a rather different model, in 
which particles aggregate, not undergo fragmentation, also studied by Takayasu earlier [25]. 


2 



FIG. 1: Structure of the fragments formed in our Model A. The picture shown is obtained from a numerical simulation, size 
of lattice L = 80, evolved to depth 150, with the crack- splitting rate A = 0.1. Note the wide distribution of fragment sizes. 


II. DEFINITION OF THE MODEL 


We start with an intact cylinder with no cracks. At the start, we introduce one or more cracks at the left boundary. 
We can define a discrete-time update rule (Model A), or continuous-time update rule ( Model B). 


Model A: We consider the model defined on a square grid, drawn on the surface of semi-infinite cylinder. The 
sites of the lattice are labelled by integers {x,t), with 1 < x < 2 L, t > 0, and (x +1 ) even. The neighbours of site 
(a;, t) are four the sites (x ± 1, t ± 1). The periodic boundary conditions along the cylinder are imposed by identifying 
(X + 2L,t) = (x,t), for all x,t. This corresponds to the axis of the cylinder being in the direction (1,1) of the lattice 

( @ • 

We consider discrete-time parallel update, where each crack moves inwards with speed 1. Then, the cracks reach 
up to distance t along the cylinder axis, and we may equivalently think of t as the time coordinate. If there is a crack 
at site {x,t), it moves according to the following rule: 

a) With probability A, it splits into two cracks, and these cracks go to (x + l,t + 1) and (x — 1 ,t + 1). Clearly, 
0 < A < 1. 

b) If the crack does not split, it moves to any one of the two forward neighbours (x + 1, t + 1) or (x — 1, t + 1) with 
equal probability. 

c) If two cracks reach a site, they merge, and move as a single crack. 


Model B: Here we consider the time-coordinate t to be a continuous variable, t > 0, but the spatial coordinate is 
still discrete, 0 < x < 2 L, with periodic boundary conditions x + 2L = x. A configuration at time t is specified by 
the occupation numbers {n^.}, where n x takes values 1 and 0, depending on whether there is a crack at x at time t or 
not. 

The configurations {n x } undergo continuous-time Markovian evolution with the following rates: 

a) A crack at position x jumps to site x' with rate 1, if x' is a nearest neighbour of x. 

b) An empty neighbour x' of an occupied site x becomes occupied with rate A, due to crack-splitting at x. 

c) If a crack arrives at a site already occupied by a crack, it merges with it, and the state of the site does not change. 
If we represent the configuration by a bit string, e.g. ...00101001110001.., the transition rates are given by the 

equations 
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FIG. 2: A system of cracks propagating into a plane leaving behind a sea of fragments in our Model A. 

(1) 
( 2 ) 


01 10, with rate 1; 

01 —>■ 11, 10 —*■ 11, with rate A; 
11 —> 01, 11 —*■ 10, with rate 1. 

Fig-! shows a schematic time-history of evolution in this case. 


III. CALCULATION OF THE STEADY STATE 

In our model, the number of cracks at any given time ( equivalently at distance along the axis of the cylinder) is 
not conserved. So long as at the start of the configuration t = 0, there is at least one crack, as the cracks propagate 
along the cylinder, they may split, or merge, and for large times, there is a non-trivial steady state of the system 
with a non-zero density p{ A) of cracks, and a finite mean area per fragment. We now determine the this steady state, 
and the asymptotic density. 


Model A: Let us label the bonds of the lattice by coordinates of their midpoints. Then a bond has coordinates 
(to + 1/2 ,n + 1/2), where to and n are any integers, with 0 < to < M — 1, and 0 < n < oo. The configuration S 
of cracks is specified by giving the occupation numbers {cr(m,n)} of all the bonds, where <r(m,n) = 1, if the bond 
(to. + l/2,n + 1/2) is occupied by a crack, and zero otherwise. The restriction of S to vertical row t will be denoted 
by E t 

We can think of the set {ct(to, n)} as the evolution history of a set of 2 M Ising spins on a line, undergoing stochastic 
Markovian dynamics. Here <t(to, n) is +1 if the m-th spin is up at time t , and zero otherwise. In this formulation, 
the Markovian evolution is defined as follows: At odd times, we consider pairs of adjacent spins [cr(2j, t),a(2j + 1, t)]. 
These evolve according to the following rules: 

a) [0,0] remains unchanged. 

b) If [ct, cr'] [0,0], it is reset to (1, 0), (0,1) or (1,1) with probabilities (1 — A)/2, (1 — A)/2, A respectively. 

At even times t, we apply the same rule, but choose the pairs to be [cr(2j + l,t),a(2j + 2, t)]. 

Now, we specify the state of the system at time f, by giving a probability measure on the space of spin configurations. 
Let Prob(E t ) denote the probability that the configuration of Ising spins at time t will be found to be S ( . Note that 
there are (2 2L — 1) allowed values of E*, as the configuration with all spins down is not reached in the steady state, 
if we start with at least one crack present at t = 0. The Master equation for the evolution of Prob(E t ) is 

Prob(E t+1 )=W[E t+1 ,E t ] Prob(E t ) (3) 

where W is the Markov transition matrix. 

We note that the transition rates here satisfy the detailed balance condition corresponding to the Hamiltonian 

M 

H(cr) = —K cr i 

i —1 


( 4 ) 
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FIG. 3: An example of the fragments of the plane left behind in model B, after the crack front has advanced beyond the 
region shown. Here time is continuous. The blue and red lines denote world lines of sites, and of cracks respectively. The green 
dots denotes crack creation, and black dots denote crack mergers. 

where e K = 2A/(1 — A). 

There are 2 2i states, and W is a 4 L x 4 L matrix. There are two steady state eigenvectors: one in which there are 
no cracks, and one in which a configuration with n cracks has a probability proportional to e Kn . 

Model B: In this case, the analysis is simpler. We note that probability of transition of the state of two adjacent 
sites from the state 01 to 11 is A, but the rate of transition from 11 to 01 is 1. Also, there is a rate 1 for the transition 
01 to 10, and vice-versa. Then, detailed balance is clearly satisfied, for the Hamiltonian function given by Eq.Q, 
with A = e K . 


IV. DISTRIBUTION OF FRAGMENT SIZES 


We first discuss the Model A. The cracks divide the plane into fragments of different sizes. We note that each 
fragment is a polygon having a unique starting point, and a unique end-point which are the boundary sites with 
minimum and maximum values of the time-coordinate respectively. The boundary itself consists of two directed 
walks. These shapes are known as staircase polygons in literature. 

The left boundary of the fragment is a simple directed biassed random walk, where the step occurs to the right 
with probability (1 + A)/2, and to the left with probability (1 — A)/2. Similarly, the right boundary is an independent 
walk with step to the left with probability (1 + A)/2, and step to the right occurs with probability (1 — A)/2. The 
fragment ends when these two boundaries meet. 

Consider the stochastic evolution of the system. Suppose there is a tip-splitting event at some site (xo,to), which 
starts a new fragment. Let Prob(A, A) be the probability that this fragment will have area A. For small A , these are 
easily determined by explicit enumeration of possible shapes. For example, consider a fragment starting at ( Xo,to )• 
We can have A = 1, only if the fragment boundary, when at the site (xo + 1, to + 1) takes the next step to the left, and 
the boundary, when at (x 0 — 1, to + 1) takes the next step to the right. Thus we have Prob(A = 1, A) = (1 + A) 2 /4. 
Similarly, Prob(A = 2, A) = (1 + A) 2 (l — A 2 )/8. More generally, if a given fragment shape S has perimeter 2P, it 
occurs with a probability (1 + A) 2 (l — X 2 ) p ~ 1 2~ p ~ 1 . 

Thus, the question of finding Prob(A) for general A is reduced to the finding the number of different staircase 
polygons of area A. This has studied earlier by Prellberg jJJO], and Prellberg and Brak [51]. Let the number of 
a staircase polygon having 2 h horizontal bonds, and 2v vertical bonds, and area A be C(h,v,A). We define the 
generating function 

G(x t y,A)='£ l * h V v 9 A ( 5 ) 

h,v,A 


Prellberg showed That G(x, y, A) satisfies a functional equation, and used it to determine G(x, y, q ) exactly. One 
gets 


G(x,y,q) = y 


' H(q 2 x, qy, q) 
H(qx,y,q) 



( 6 ) 
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where 


with 


H{x,y,q) 


E 


n —0 


(-x) n q&) 

(q',q)n{y;q)n 


(r;q) n = ^ _ tqm ) 

m =0 


(7) 


( 8 ) 


These explicit formulas are a bit complicated. These simplify in the continuum limit, corresponding to A —> 0. 
The continuum problem of determining the probability distribution of area in a Brownian bridge has been studied in 
[S3 155] - They showed that in this case, the distribution function Prob(A,A) is the Airy distribution. For small A, 
Prob(A, A) has the scaling form 

Prob(A,A) ~ A 4 ff (A 3 A). (9) 

In the limit A —> 0, and fixed A, Prob(A, A) tends to a non-zero value. In this case, the boundaries of the fragment 
do unbiased random walk,and the corresponding probabilities are exactly as found in the directed abelian sandpile 
model in 2-dinrensions m- For this case, simple dimensional arguments already show that Prob(A) ~ A 4 / 3 37 ] • 
This implies that 


Prob(A, A —» 0) = Prob(A) J 4SM (10) 

where Prob(A)^gM is the probability that adding a grain in the steady state of 2-dimensional directed abelian sandpile 
model will cause an avalanche with exactly A topplings E3HS |. It is known that Prob(A) y 4gM ~ A 4 / 3 for large A, 
we must have 


g(x) ~ Kx 4//3 , for small x. (11) 

For large x, the Airy distribution, g(x) decreases as exp(—a; 1 / 2 ) [33} . In the literature, there has been some discussion 
about the size-distribution having exponential or stretched exponential decay for large sizes [I] . 

Note the unusual A 4 scaling of the prefactor in Eq. ([9]). This arises because the function g(x), which varies as x ~ 4 / 3 
for small x is not integrable. Then, if we take the scaling limit 


^Prob(A,A) 

A =1 



dxg(x), 


( 12 ) 


we have to impose a cut-off of x m i n = 0( A 3 ) on the lower limit of the integral. Then the integral diverges as x^ nln , 
which cancels the extra power of A in the Eq. ([9|. 

For finite but small A, the mean density of cracks varies as 1/A, and the average longitudinal extent of a fragment 
scales as 1/A 2 . Thus, the mean size of a fragment varies as 1/A 3 . 

Now, we consider Model B. Let Prob(A, m)dA be the probability that a segment of width to will generate a fragment 
with additional area between A and A + dA. The Laplace transform of this distribution will be denoted by P(s,m), 
and is defined by 

/»oo 

P(s,m) = dA Prob(A, m)e~ sA (13) 

Jo 

Consider the time evolution of a fragment whose transverse size is to at some time t = to. In time, the boundaries of 
this fragment will evolve stochastically, till they merge. The segment of length to can become segment of length mil, if 
either end makes a diffusive jump ( with rate 2 for each), or if either end creates a new crack which decreases the size of 
this segment (total rate 2A). The waiting time t\ before any one of these occurs is an exponentially distributed random 
variable. The probability that this occurs between time t\ and t\ + dt\ given by Prob(fi)dti = (4 + 2A)e“( 4+2A ) tl <i<i. 
The total area of the fragment is sum of two independent random variables: the area A± = mt\ up to this time, 
and remainder area A 2 . The remainder has starting length to + 1 with probability p + = 1/(2 + A), and to — 1, with 
probability p_ = (1 + A)/(2 + A). This gives 


P(a, to ) 


Prob(ti)dt\e 


—smt 1 


U 0 


p + P(s,m + 1) + p_ P(s, 1 


1) 


(14) 
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which simplifies to 

(4 + 2A + ms)P(s, to) = 2 P(s, m + 1) + (2 + 2A )P(s, m — 1) (15) 

Write P(sm) = (1 + A )~ s ^ 2 Q(s , to). Then Q(s,m) satisfies the equation of the form 
4 “I - 2A “I - 777 S 

(- j — 2 )Q(s, to) = [Q(s, to + 1) + Q(s, m - 1) - 2 Q(s, to)] (16) 

v 1 + A 

which is a discrete variant of the Airy equation. In the continuum limit of small A , it reduces to the Airy equation, 
and the scaling form of the fragment distribution is the same as studied in [32] . 


V. ULTRA-LOCALITY OF CORRELATIONS BETWEEN FRAGMENTS 

This model has the remarkable property that two fragment starting at any two sites (x\,t\) and [x 2 ,t- 2 ) are 
completely uncorrelated, unless they are actually touching each other. This is true not only for their area, but also 
their shapes. This lack of correlations between fragments is very surprizing, in view of the fact that the condition 
that the arrangement of fragments has to tile the plane imposes very strong geometrical constraints on their shapes. 
The probability that a fragment starting at (x±, tp) will have a given have any given shape S, depends only on S, and 
is completely independent of the state of other sites at time t\. If there is another crack that comes near this crack, 
at best it can merge with this fragment’s boundary, but this will not influence its subsequent motion. 

We have hitherto considered all cracks moving into the unfragmented region with unit speed. This assumption is 
not really necessary. At the start of the simulation, we can assign to each site a random variable, taking one of three 
possible possible values, which tells what would happen if a crack reaches that site: will it split into two, or go left 
unsplit, or go right unsplit. Then, we can start with any configuration of cracks at the boundary on the left, and 
evolve cracks in any order, and the final configuration is independent of the order in which the cracks are updated, 
and the waiting time between these updates. 

This kind of “ultra-locality” of cluster shapes is also seen in percolation models, where the shapes of two percolation 
clusters are perfectly uncorrelated, unless they have a common boundary. The present model may also be seen as 
a model of this type. This property is much less obvious, if we look at the transfer matrix of the model. Consider 
the configuration of a particular column t = to- We assume that a crack-splitting occurs at a point (xo,to)- The 
probability that the generated fargment with (xo,to) as its left-most point has area A, denoted by Prob(A, A) above, 
would in general be expressible in terms the eigenvalues and eigenvectors of the transfer matrix, and matrix elements 
of some projection operators between eigenvectors. That this probability does not depend on the starting vector, so 
long as a splitting occurs at (xo,to), requires the transfer matrix to have very special structure. For example, the 
transfer matrix W can be written as 


W = TiT 2 (17) 

where Ti and T 2 are the matrices cooresponding to evolution at odd and even times. But T 2 can be written as a 
direct product of L 4 x 4 matrices, where each matrix gives the transition probabilities from the state of a pair of 
adjacent sites cr(2j — 1, t), cr(‘2j, t) to the state a(2j — 1, t + l),a(2j, t + 1). Each of these is the configuration basis is 
easily seen to be of the form 


T 0 0 0' 

0 a a a 
0 a a a 
. 0666 . 

where a = (1 — A)/2, 6 = A. Clearly, the rank of T 2 , and hence of W is at most 2 L . Out of the A L eigenvalues 
of L , at least 4 L — 2 L are exactly zero. The large number of zero eigenvalues provides the mathematical mechanism 
underlying the the ultra-locality in this model. 

I thank S. N. Majumdar for discussions, and for pointing out ref. m to me. This work is supported in part by the 
Department of Science and Technology, Government of India via grant DST-SR/S2/JCB-24/2005. 
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